↜ Back to index Introduction to Numerical Analysis 2

Part b–Lecture 1

Goal of this quarter

During this quarter we will learn methods to solve Ax = b that are faster than the Gaussian elimination (and even the LU decomposition) for special matrices. The main goal is to implement the conjugate gradient method (共役勾配法) for symmetric positive definite matrices. The basic code will turn out to be easier than the Gaussian elimination.

But first we learn how to define and use functions in Fortran. It will be more convenient to write the code using functions.

Functions in Fortran

In mathematics, a function (関数) like \sin or \operatorname{abs} is a rule x \mapsto \sin x that given one or more inputs produces an output. We have to specify the domain and the range, that is, what type of inputs are allowed and what type of output is to be expected: \sin: \mathbb R \to \mathbb R.

A function in Fortran captures a similar idea. In the function definition, we specify the name of the function, the types of inputs and the output, and provide the code that produces the output based on the inputs. The inputs are also called function arguments. We then can use the function in other parts of the code to easily get the output for the given inputs.

The main advantage of using functions is that we do not need to remember how exactly the function is implemented to use it.

To illustrate the syntax, let us define a function add_two that simply adds 2 to its argument and returns the new value, and a function foo that multiplies its two arguments increased by two:

implicit none

print *, add_two(2.0)
print *, add_two(10.0)
print *, foo(1.0, 3.0)

contains
    function add_two(x)
        real x
        real add_two

        add_two = x + 2
    end

    function foo(x, y)
        real x, y, foo

        foo = add_two(x) * add_two(y)
    end
end

Before the end command that concludes the program we included a new section contains where one can define functions that are then accessible from the main program (the part before contains).

A function definition starts with the keyword function and ends with the keyword end. function must be followed by the name of the function and the list of arguments in parentheses, separated by commas. The lines that follow then declare all the variables used in the function code, including the types of the arguments (x in the above code) and the output (add_two, the same as the name of the function). The rest of the code before the end of the function is the code used to produce the output value.

We can call functions from other functions.

Here is a more complicated example of a function that computes a factorial.

implicit none
integer i

do i = 1, 5
    print *, 'Factorial of', i, 'is', fact(i)
end do

contains
    function fact(n)
        integer n, fact, i

        fact = 1
        do i = 1, n
            fact = fact * i
        end do
    end
end

This time we are using a local variable i within the function. Note one very important fact: since i is declared inside the function fact, it is a different variable from the variable i in the main program. In particular, the do loop inside fact will not change the value of i in the main program and the produced result is correct.

Note that the return variable fact can be used like any other variable. The value that is returned is the value of the variable when the function ends.

Exercise 1. Try removing i from the line integer n, fact, i and rerun the program.

Exercise 2. Implement a function maximum that takes two real numbers and returns their maximum. Implement another function positive_part that takes one real number x as an input and returns \max(x, 0) using the function maximum.

Using arrays as function arguments

Functions in Fortran can accept arrays as arguments and even return arrays.

The following function takes a one dimensional array of any size (indicated by (:) in the declaration of x) and returns the array of the differences of the consecutive elements.

implicit none
real a(3), b(5)

a = [3.0, 2.5, 5.0]
b = [1.0, 3.0, 2.0, 2.5, 5.0]

print *, diffs(a)
print *, diffs(b)

contains
    function diffs(x)
        real x(:)
        real diffs(size(x) - 1)
        integer i

        do i = 1, size(diffs)
            diffs(i) = x(i+1) - x(i)
        end do  
    end
end

Here is a more complicated function that takes a matrix as an input and returns the array of the sums of each row of the matrix.

implicit none
real a(3,2)

! initialize a by setting its rows
a(1, :) = [1., -1.]
a(2, :) = [2., 3.]
a(3, :) = [6., 7.]

print *, row_sum(a)

contains
    function row_sum(a)
        real a(:,:)                 ! matrix of any size
        real row_sum(size(a, 1))    ! size equal to the number of rows of a
        integer i, j

        ! iterate of rows
        do i = 1, size(a, 1)
            row_sum(i) = 0
            ! iterate over columns
            do j = 1, size(a, 2)
                row_sum(i) = row_sum(i) + a(i, j)
            end do
        end do
    end
end

The code takes as input a two dimension array of any size indicated as (:,:), and returns a one-dimensional array with the length equal to the number of rows in a.

Built-in function `size`

You can read more about the built-in function size here. In short, size(a) returns the number of all the elements of array a, while size(a, 1) returns the number of rows and size(a, 2) the number of columns in the matrix represented by a.

real a(5, 3)

print *, size(a, 1)
print *, size(a, 2)
print *, size(a)

end

prints

           5
           3
          15

Exercise 3.

Implement a function inner that takes two one-dimensional arrays x and y of the same size representing two vectors x, y and returns their inner product x \cdot y = (x, y) := x_1 y_1 + \cdots + x_n y_n.

Given arrays

real a(3), b(3)

a = [3.0, 2.5, 5.0]
b = [1.0, -3.0, 2.0]

print *, inner(a, b)

we get

   5.50000000

Note. Fortran 95 already has a built-in function dot_product that does the same thing.

Exercise 4.

Implement a function mv_mul that takes a two-dimensional array a representing a matrix A and a one-dimensional array x representing a column vector x as input and return the one-dimensional array that contains the component of the product Ax.

Hint. Recall that (Ax)_i = a_{i1} x_1 + a_{i2} x_2 + \cdots + a_{in} x_n. This is the inner product of row i of matrix A and vector x. You can therefore use the function inner from the previous exercise if you want.

With the following arrays:

real a(3,2), x(2)

a(1, :) = [1., -1.]
a(2, :) = [2., 3.]
a(3, :) = [6., 7.]

x = [3., -1.]

print *, mv_mul(a, x)

we get

   4.00000000       3.00000000       11.0000000

Exercise 5.

Implement a function norm that that takes a one-dimensional array x and returns the Euclidean norm \|x\| := \sqrt{x \cdot x} of the represented vector x.

Hint. You can use the function inner implemented before.

For array

real x(3)

x = [1., 2., -3.]

print *, norm(x)

we get

   3.74165750

Note. Fortran 2008 has a built-in function norm2 that does the same thing.

More on Fortran functions

Recursive functions

In programming, functions are often called recursively (再帰呼び出し): a function calls itself with different argument values.

An important example is the definition of a factorial (階乗): n! = \begin{cases} 1 & n = 0,\\ n \cdot (n-1)! & n > 0. \end{cases}

In Fortran, we had to declare the function as recursive. In most other languages, all functions can be called recursively. Fortran is very old…

implicit none

print *, fact(5)
print *, fact(12)

contains
 ! 👇👇👇👇👇
    recursive function fact(n) result(r)
        integer n, r

        if (n == 0) then
            r = 1
            return              ! return from the function
        endif 

        r = n * fact(n - 1)
    end
end

We also had to specify the name for a return value using result(r).

return command is useful if you want to return from a function early.

Exercise 6. Write a function fib(n) for computing the Fibonacci number F_n using the recursive formula F_n := \begin{cases} 0, & n = 0,\\ 1, & n = 1,\\ F_{n-2} + F_{n-1}, & n \geq 2. \end{cases}

How many times is fib called when computing fib(5)? How many times for fib(10)?

Exercise 7. Write a function gcd(n, k) that computes the greatest common divisor (最大公約数) of two positive integers n and k recursively using the Euclidean algorithm (ユークリッドの互除法):

gcd(n, 0) is n
if k > 0, then gcd(n, k) is the same as gcd(k, mod(n, k)).

mod(n, k) is the Fortran’s built-in function for the remainder (剰余).

Extra exercises

Exercise 8. Write a program that attempts to find the eigenvector of a given symmetric matrix A corresponding to the largest eigenvalue using the following algorithm.

Set x^{(0)} to some nonzero vector.
Iterate for i = 1, 2, …
- Compute x^{(i)} = \frac{Ax^{(i-1)}}{\|Ax^{(i-1)}\|} Here the program should report an error if \|Ax^{(i-1)}\| is too small, say 10^{-6}.
- Print the value of x^{(i)} to track progress.
- Stop the iteration if \|x^{(i)} - x^{(i-1)}\| < \varepsilon.

Use \varepsilon = 10^{-6}.

Hints.

Use the functions mv_mul and norm that you implemented in the previous exercises.
In Fortran, you can divide all the elements of an array x by a real number y by simply writing x / y. No need to use do loops. Similarly, to compute the element-wise difference of two arrays v and w of the same size you can use v - w.
In numerical algorithms we usually restrict the maximum number of iterations in case the algorithm does not converge. You can use for example do i = 1, 100 to limit the number of iterations to 100.

Sample run of the program with

real a(3,3), x(3)

a(1,:) = [1.0, 2.0, 3.0]
a(2,:) = [2.0, 1.0, 4.0]
a(3,:) = [3.0, 4.0, 1.0]

x = [1., 0., 0.]

  0.267261237      0.534522474      0.801783681    
  0.549971938      0.628539383      0.549971938    
  0.490831792      0.557763398      0.669316173    
  0.511315823      0.596535206      0.618629038    
  0.503328860      0.578738987      0.641655087    
  0.506864071      0.586905181      0.631372452    
  0.505297840      0.583228707      0.636017621    
  0.506004035      0.584889233      0.633927822    
  0.505686522      0.584141433      0.634870052    
  0.505829632      0.584478498      0.634445608    
  0.505765200      0.584326684      0.634636879    
  0.505794227      0.584395111      0.634550691    
  0.505781174      0.584364235      0.634589553    
  0.505787075      0.584378183      0.634572029    
  0.505784392      0.584371865      0.634579957    
  0.505785584      0.584374666      0.634576380    
  0.505785048      0.584373474      0.634577930    
  0.505785286      0.584374011      0.634577274

Goal of this quarter

Functions in Fortran

Using arrays as function arguments

Built-in function size

More on Fortran functions

Recursive functions

Extra exercises

Built-in function `size`